stanEx5
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]'
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}hierarchical model
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
scale_shape_manual(values=1:k)
qplot(x,y,col=c)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionestimate as no class
ex8-0.stan
y~N(b00+b10*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)## Initial log joint probability = -4279.57
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -147.802 0.00185173 0.00427274 0.9623 0.9623 34
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -147.80
## b0 8.30
## b1 3.06
## s 11.66
## m0[1] 12.68
## m0[2] 54.80
## m0[3] 50.94
## m0[4] 51.73
## m0[5] 40.74
## m0[6] 58.71
## m0[7] 14.48
## m0[8] 21.19
## m0[9] 42.97
## m0[10] 57.17
## m0[11] 27.00
## m0[12] 61.03
## m0[13] 60.31
## m0[14] 32.99
## m0[15] 19.88
## m0[16] 51.52
## m0[17] 26.94
## m0[18] 23.45
## m0[19] 67.97
## m0[20] 57.35
## m0[21] 15.96
## m0[22] 32.89
## m0[23] 18.16
## m0[24] 65.77
## m0[25] 16.15
## m0[26] 37.52
## m0[27] 14.73
## m0[28] 65.51
## m0[29] 47.44
## m0[30] 13.93
## m0[31] 25.28
## m0[32] 36.12
## m0[33] 47.17
## m0[34] 60.20
## m0[35] 17.08
## m0[36] 36.67
## m0[37] 39.22
## m0[38] 28.08
## m0[39] 8.73
## m0[40] 24.25
## m0[41] 38.03
## m0[42] 32.88
## m0[43] 69.21
## m0[44] 20.78
## m0[45] 18.62
## m0[46] 52.52
## m0[47] 64.15
## m0[48] 33.00
## m0[49] 54.09
## m0[50] 39.93
## y0[1] 14.25
## y0[2] 61.30
## y0[3] 54.72
## y0[4] 61.84
## y0[5] 48.18
## y0[6] 67.23
## y0[7] 8.33
## y0[8] 12.94
## y0[9] 35.61
## y0[10] 40.00
## y0[11] 22.85
## y0[12] 34.27
## y0[13] 65.45
## y0[14] 20.48
## y0[15] 43.05
## y0[16] 51.08
## y0[17] 27.82
## y0[18] 14.98
## y0[19] 58.33
## y0[20] 59.86
## y0[21] 5.64
## y0[22] 34.98
## y0[23] 14.57
## y0[24] 64.19
## y0[25] 3.89
## y0[26] 28.39
## y0[27] 29.21
## y0[28] 68.45
## y0[29] 40.53
## y0[30] 16.80
## y0[31] 27.90
## y0[32] 46.10
## y0[33] 41.44
## y0[34] 62.64
## y0[35] 23.20
## y0[36] 32.50
## y0[37] 36.19
## y0[38] 42.13
## y0[39] 7.30
## y0[40] 11.17
## y0[41] 43.63
## y0[42] 30.78
## y0[43] 64.93
## y0[44] 30.09
## y0[45] 20.63
## y0[46] 87.02
## y0[47] 66.83
## y0[48] 52.52
## y0[49] 21.48
## y0[50] 41.61
## m1[1] 8.73
## m1[2] 15.45
## m1[3] 22.17
## m1[4] 28.89
## m1[5] 35.61
## m1[6] 42.33
## m1[7] 49.05
## m1[8] 55.77
## m1[9] 62.49
## m1[10] 69.21
## y1[1] -18.31
## y1[2] 31.64
## y1[3] 29.22
## y1[4] 19.51
## y1[5] 44.08
## y1[6] 42.64
## y1[7] 55.84
## y1[8] 53.35
## y1[9] 59.47
## y1[10] 74.10
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -146.88 -146.57 1.24 1.00 -149.47 -145.54 1.00 631 755
## b0 8.45 8.50 3.46 3.53 2.80 13.90 1.01 279 265
## b1 3.04 3.04 0.30 0.31 2.56 3.54 1.01 355 507
## s 12.21 12.08 1.28 1.22 10.34 14.44 1.00 1923 1288
## m0[1] 12.80 12.83 3.09 3.17 7.71 17.71 1.01 280 274
## m0[2] 54.73 54.79 2.36 2.33 50.81 58.51 1.00 1351 974
## m0[3] 50.89 50.93 2.11 2.05 47.34 54.25 1.00 1712 1126
## m0[4] 51.68 51.73 2.16 2.09 48.02 55.16 1.00 1699 1086
## m0[5] 40.74 40.71 1.73 1.68 37.82 43.49 1.00 930 1163
## m0[6] 58.62 58.62 2.64 2.65 54.23 62.80 1.00 1068 928
## m0[7] 14.60 14.64 2.94 3.01 9.82 19.29 1.01 281 288
## m0[8] 21.28 21.25 2.43 2.47 17.33 25.16 1.01 293 355
## m0[9] 42.96 42.97 1.77 1.71 40.02 45.79 1.00 1217 1198
## m0[10] 57.09 57.12 2.52 2.52 52.92 61.10 1.00 1165 986
## m0[11] 27.06 27.05 2.06 2.09 23.65 30.40 1.01 327 418
## m0[12] 60.93 60.92 2.82 2.82 56.22 65.44 1.00 948 938
## m0[13] 60.21 60.20 2.76 2.75 55.59 64.62 1.00 981 928
## m0[14] 33.02 33.01 1.80 1.77 30.02 35.93 1.01 424 641
## m0[15] 19.98 19.98 2.52 2.57 15.86 23.99 1.01 289 349
## m0[16] 51.47 51.52 2.15 2.08 47.84 54.92 1.00 1703 1099
## m0[17] 27.00 27.00 2.06 2.09 23.59 30.34 1.01 326 418
## m0[18] 23.53 23.51 2.27 2.29 19.79 27.17 1.01 302 408
## m0[19] 67.85 67.87 3.39 3.37 62.14 73.37 1.00 739 876
## m0[20] 57.27 57.31 2.54 2.52 53.10 61.29 1.00 1153 986
## m0[21] 16.07 16.12 2.82 2.92 11.48 20.58 1.01 282 313
## m0[22] 32.92 32.91 1.80 1.78 29.91 35.83 1.01 421 626
## m0[23] 18.26 18.28 2.65 2.72 13.95 22.46 1.01 285 334
## m0[24] 65.65 65.69 3.21 3.20 60.32 70.78 1.00 791 874
## m0[25] 16.26 16.30 2.81 2.90 11.69 20.75 1.01 283 313
## m0[26] 37.53 37.51 1.71 1.67 34.66 40.30 1.01 638 970
## m0[27] 14.85 14.89 2.92 3.00 10.11 19.52 1.01 282 288
## m0[28] 65.40 65.44 3.18 3.17 60.11 70.48 1.00 797 874
## m0[29] 47.41 47.45 1.93 1.90 44.10 50.49 1.00 1720 1265
## m0[30] 14.05 14.09 2.98 3.05 9.17 18.80 1.01 281 281
## m0[31] 25.34 25.34 2.16 2.18 21.79 28.83 1.01 312 421
## m0[32] 36.14 36.10 1.73 1.70 33.24 38.92 1.01 545 930
## m0[33] 47.14 47.17 1.92 1.89 43.86 50.19 1.00 1715 1287
## m0[34] 60.10 60.09 2.75 2.75 55.49 64.50 1.00 985 928
## m0[35] 17.18 17.22 2.73 2.80 12.69 21.52 1.01 284 311
## m0[36] 36.68 36.64 1.72 1.69 33.76 39.46 1.01 578 940
## m0[37] 39.22 39.20 1.71 1.67 36.32 41.97 1.00 782 1027
## m0[38] 28.14 28.16 2.00 2.04 24.77 31.40 1.01 336 482
## m0[39] 8.88 8.94 3.42 3.49 3.29 14.27 1.01 279 266
## m0[40] 24.32 24.32 2.22 2.25 20.68 27.89 1.01 306 412
## m0[41] 38.04 38.01 1.71 1.69 35.13 40.80 1.00 680 994
## m0[42] 32.92 32.91 1.80 1.78 29.91 35.82 1.01 421 626
## m0[43] 69.08 69.10 3.50 3.46 63.20 74.77 1.00 716 866
## m0[44] 20.87 20.84 2.46 2.50 16.85 24.78 1.01 291 366
## m0[45] 18.72 18.73 2.62 2.67 14.47 22.88 1.01 286 336
## m0[46] 52.46 52.53 2.21 2.16 48.75 56.02 1.00 1678 1016
## m0[47] 64.03 64.08 3.07 3.07 58.95 68.97 1.00 838 908
## m0[48] 33.03 33.02 1.80 1.77 30.03 35.94 1.01 424 641
## m0[49] 54.02 54.10 2.31 2.26 50.16 57.73 1.00 1459 947
## m0[50] 39.93 39.92 1.72 1.67 37.03 42.68 1.00 847 1141
## y0[1] 13.35 13.54 12.52 12.37 -7.57 33.51 1.00 1449 1805
## y0[2] 55.12 55.08 12.38 12.19 35.24 75.10 1.00 2165 2105
## y0[3] 50.92 50.63 12.63 12.04 30.66 72.00 1.00 1995 1856
## y0[4] 51.53 51.47 12.65 12.50 31.08 72.46 1.00 1869 1804
## y0[5] 40.84 40.52 12.56 12.34 20.60 61.41 1.00 2014 1931
## y0[6] 58.49 58.73 12.22 11.92 38.19 78.12 1.00 2011 1906
## y0[7] 14.26 14.56 12.60 12.81 -6.51 35.05 1.00 1548 1931
## y0[8] 21.31 21.15 12.53 11.96 0.18 41.57 1.00 1797 1972
## y0[9] 43.31 43.28 12.43 12.26 23.21 63.93 1.00 1929 1798
## y0[10] 56.95 57.38 12.22 12.20 36.23 76.93 1.00 1911 1973
## y0[11] 26.92 26.68 12.37 12.01 6.66 47.09 1.00 2004 1890
## y0[12] 60.91 60.98 12.25 11.90 40.97 80.92 1.00 1843 1900
## y0[13] 60.49 60.43 12.49 12.34 40.02 81.11 1.00 1990 2058
## y0[14] 33.42 33.81 12.22 12.14 12.51 53.44 1.00 1959 1996
## y0[15] 20.25 20.89 12.50 12.76 -0.67 40.87 1.00 1998 1828
## y0[16] 51.24 51.16 12.41 12.30 31.41 71.70 1.00 1650 2008
## y0[17] 27.07 27.03 12.21 12.44 6.68 46.84 1.00 2015 1933
## y0[18] 23.47 23.62 12.50 12.26 2.88 43.42 1.00 2011 1934
## y0[19] 67.53 67.60 12.95 12.95 46.69 88.10 1.00 1919 1972
## y0[20] 57.14 57.31 12.80 12.66 36.19 77.47 1.00 1773 1969
## y0[21] 15.67 16.04 12.69 12.35 -5.39 36.31 1.00 1587 1757
## y0[22] 33.13 32.94 12.31 12.14 12.82 53.35 1.00 2131 2082
## y0[23] 18.27 18.24 12.13 11.75 -1.72 38.37 1.00 1910 1779
## y0[24] 66.09 65.73 12.83 12.57 44.54 87.55 1.00 1810 1795
## y0[25] 15.95 15.94 12.26 11.73 -4.13 36.33 1.00 1679 1918
## y0[26] 37.10 36.98 12.74 12.05 15.71 58.49 1.00 1753 1882
## y0[27] 15.13 14.93 12.64 12.41 -5.68 35.15 1.00 1733 1741
## y0[28] 65.16 65.09 12.80 12.13 43.79 86.44 1.00 1744 1910
## y0[29] 47.08 46.91 12.97 12.94 25.67 68.26 1.00 2042 1779
## y0[30] 13.75 14.07 12.73 12.92 -7.48 34.61 1.00 1816 1766
## y0[31] 25.29 25.18 12.70 12.51 4.56 46.49 1.00 1719 1709
## y0[32] 36.55 36.79 12.61 12.82 15.45 57.82 1.00 1927 2057
## y0[33] 47.02 47.34 12.52 12.92 26.37 67.02 1.00 2116 1882
## y0[34] 60.43 60.45 12.58 12.39 39.31 80.96 1.00 1845 1907
## y0[35] 17.17 16.99 12.48 12.15 -3.18 37.83 1.00 1899 1929
## y0[36] 36.49 36.58 12.06 12.12 16.42 56.75 1.00 1936 1659
## y0[37] 39.27 39.16 12.21 11.99 19.05 58.70 1.00 2013 1818
## y0[38] 28.09 28.56 12.18 12.22 7.97 47.40 1.00 1911 1923
## y0[39] 9.07 8.89 12.96 12.50 -11.68 30.28 1.00 1314 1725
## y0[40] 23.97 23.87 12.62 12.00 3.71 45.07 1.00 1743 1719
## y0[41] 38.17 38.15 12.19 12.25 18.53 58.39 1.00 1810 1812
## y0[42] 33.04 32.97 12.21 12.35 13.06 53.00 1.00 1815 1794
## y0[43] 68.92 68.78 13.23 13.36 47.46 90.55 1.00 1976 1884
## y0[44] 20.60 20.62 12.49 12.24 0.27 41.04 1.00 1819 1932
## y0[45] 18.07 18.07 12.56 12.48 -2.96 38.29 1.00 1934 1861
## y0[46] 52.65 52.69 12.26 12.83 33.12 72.39 1.00 1888 1971
## y0[47] 64.00 63.88 12.82 12.63 43.21 85.26 1.00 1972 1946
## y0[48] 32.88 32.61 12.08 11.71 12.92 52.86 1.00 2129 1933
## y0[49] 53.50 53.88 12.39 12.35 32.66 72.97 1.00 1986 2014
## y0[50] 40.28 40.31 12.34 12.43 20.42 60.31 1.00 1705 1838
## m1[1] 8.88 8.94 3.42 3.49 3.29 14.27 1.01 279 266
## m1[2] 15.57 15.63 2.86 2.95 10.92 20.15 1.01 282 307
## m1[3] 22.26 22.23 2.36 2.39 18.40 26.03 1.01 296 369
## m1[4] 28.94 28.96 1.96 1.99 25.65 32.16 1.01 345 528
## m1[5] 35.63 35.59 1.74 1.70 32.72 38.44 1.01 517 857
## m1[6] 42.32 42.32 1.75 1.70 39.39 45.15 1.00 1126 1224
## m1[7] 49.01 49.05 2.01 1.96 45.58 52.20 1.00 1725 1196
## m1[8] 55.70 55.74 2.42 2.40 51.67 59.58 1.00 1264 974
## m1[9] 62.39 62.42 2.93 2.93 57.49 67.11 1.00 890 920
## m1[10] 69.08 69.10 3.50 3.46 63.20 74.77 1.00 716 866
## y1[1] 8.82 9.13 12.72 12.63 -11.69 29.96 1.00 1344 1656
## y1[2] 15.81 15.99 12.51 11.84 -5.60 36.09 1.00 1735 1715
## y1[3] 22.23 21.97 12.49 12.22 1.96 43.32 1.00 1841 1683
## y1[4] 28.70 28.78 12.33 12.30 8.18 48.76 1.00 1565 2015
## y1[5] 36.00 35.29 12.12 11.85 17.15 56.29 1.00 1625 1792
## y1[6] 42.29 42.36 12.37 12.33 22.01 61.89 1.00 1852 1927
## y1[7] 49.33 49.74 12.21 12.26 28.18 68.34 1.00 1892 1969
## y1[8] 55.34 54.93 12.65 12.52 34.41 76.11 1.00 1942 1740
## y1[9] 62.55 62.50 12.80 12.68 41.69 83.36 1.00 2011 1860
## y1[10] 68.85 69.12 12.68 12.62 48.24 89.15 1.00 1834 1786
estimate as independent class
all b0l,b1l do not have a distribution
class c=1~k ex12-1.stan
yc~N(b0c+b1c*x,s)
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
vector[K] b0;
vector[K] b1;
real<lower=0> s;
}
model{
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-1.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -587644
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 99 -86.3037 0.232138 0.816477 1 1 154
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 -85.7831 0.0876634 0.790418 0.2119 1 265
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 299 -85.4373 0.00807286 0.0674109 0.9062 0.9062 370
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 399 -85.4126 0.00123051 0.00769343 1 1 481
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 499 -85.4118 0.000725348 0.00279152 1 1 587
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -85.41
## b0[1] -2.33
## b0[2] -0.42
## b0[3] 26.04
## b0[4] 5.46
## b0[5] 16.22
## b0[6] 1.71
## b0[7] 12.72
## b0[8] 11.71
## b0[9] 13.13
## b1[1] 0.87
## b1[2] 4.00
## b1[3] 3.55
## b1[4] 3.99
## b1[5] 3.19
## b1[6] 3.96
## b1[7] 2.55
## b1[8] 2.21
## b1[9] 2.89
## s 3.35
## m0[1] 17.27
## m0[2] 51.53
## m0[3] 53.47
## m0[4] 43.07
## m0[5] 6.88
## m0[6] 11.98
## m0[7] 22.66
## m0[8] 16.45
## m0[9] 36.75
## m0[10] 67.15
## m0[11] 24.04
## m0[12] 74.32
## m0[13] 70.42
## m0[14] 33.66
## m0[15] 24.09
## m0[16] 76.26
## m0[17] 29.80
## m0[18] 25.36
## m0[19] 69.58
## m0[20] 53.65
## m0[21] 11.62
## m0[22] 36.39
## m0[23] 22.46
## m0[24] 67.50
## m0[25] 15.71
## m0[26] 37.81
## m0[27] 18.08
## m0[28] 53.02
## m0[29] 39.97
## m0[30] -0.73
## m0[31] 29.19
## m0[32] 35.98
## m0[33] 50.44
## m0[34] 49.18
## m0[35] 20.04
## m0[36] 39.97
## m0[37] 40.03
## m0[38] 27.31
## m0[39] 13.54
## m0[40] 26.03
## m0[41] 37.53
## m0[42] 29.46
## m0[43] 70.75
## m0[44] 17.86
## m0[45] 18.94
## m0[46] 62.30
## m0[47] 90.92
## m0[48] 33.33
## m0[49] 63.94
## m0[50] 40.97
## y0[1] 15.37
## y0[2] 44.33
## y0[3] 54.84
## y0[4] 46.67
## y0[5] -0.74
## y0[6] 14.73
## y0[7] 25.49
## y0[8] 14.23
## y0[9] 28.75
## y0[10] 72.09
## y0[11] 21.59
## y0[12] 76.05
## y0[13] 69.59
## y0[14] 31.42
## y0[15] 28.80
## y0[16] 77.22
## y0[17] 29.13
## y0[18] 31.56
## y0[19] 67.90
## y0[20] 56.16
## y0[21] 10.52
## y0[22] 34.49
## y0[23] 22.78
## y0[24] 70.25
## y0[25] 13.85
## y0[26] 37.22
## y0[27] 22.34
## y0[28] 54.94
## y0[29] 41.55
## y0[30] 4.06
## y0[31] 28.63
## y0[32] 35.30
## y0[33] 45.37
## y0[34] 45.88
## y0[35] 24.26
## y0[36] 40.22
## y0[37] 42.55
## y0[38] 24.25
## y0[39] 15.27
## y0[40] 26.94
## y0[41] 36.98
## y0[42] 24.45
## y0[43] 70.04
## y0[44] 16.32
## y0[45] 12.37
## y0[46] 60.08
## y0[47] 89.04
## y0[48] 36.06
## y0[49] 60.55
## y0[50] 34.71mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 1.7 seconds.
## Chain 3 finished in 1.8 seconds.
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 1.9 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 2.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.8 seconds.
## Total execution time: 2.1 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -96.02 -95.58 3.74 3.79 -102.85 -90.47 1.01 580 988
## b0[1] -2.24 -2.18 4.61 4.52 -9.67 5.26 1.00 2567 1463
## b0[2] -0.22 -0.21 5.99 5.84 -10.06 9.45 1.01 1968 1154
## b0[3] 27.80 26.24 24.69 24.65 -10.91 70.14 1.01 271 363
## b0[4] 5.51 5.55 3.38 3.20 -0.03 10.99 1.01 3016 1517
## b0[5] 16.23 16.19 5.04 5.12 7.81 24.51 1.01 2161 1363
## b0[6] 1.82 1.83 5.59 5.61 -7.43 11.23 1.00 1842 1302
## b0[7] 12.71 12.71 2.94 2.87 7.82 17.55 1.00 2421 1149
## b0[8] 11.43 11.63 7.24 7.21 -0.76 22.69 1.00 1331 1016
## b0[9] 13.18 13.23 2.08 2.09 9.79 16.61 1.00 2389 1199
## b1[1] 0.85 0.85 0.40 0.40 0.18 1.49 1.00 2459 1664
## b1[2] 3.99 3.98 0.64 0.62 2.96 5.06 1.01 1986 1096
## b1[3] 3.44 3.54 1.52 1.56 0.84 5.83 1.01 275 416
## b1[4] 3.99 3.99 0.37 0.36 3.40 4.60 1.00 2744 1419
## b1[5] 3.19 3.19 0.36 0.36 2.60 3.80 1.01 2145 1490
## b1[6] 3.95 3.93 0.97 0.94 2.36 5.56 1.00 1804 1388
## b1[7] 2.56 2.56 0.31 0.30 2.04 3.07 1.00 2506 1344
## b1[8] 2.23 2.22 0.51 0.52 1.44 3.10 1.00 1360 1081
## b1[9] 2.89 2.88 0.18 0.18 2.60 3.18 1.00 2231 1580
## s 4.33 4.27 0.56 0.55 3.53 5.35 1.00 1182 1456
## m0[1] 17.31 17.33 1.89 1.88 14.22 20.43 1.00 2370 1203
## m0[2] 51.61 51.64 2.73 2.59 46.99 56.14 1.00 2189 1847
## m0[3] 53.46 53.45 1.53 1.49 51.00 55.91 1.00 2175 1338
## m0[4] 43.15 43.18 1.80 1.75 40.09 46.04 1.01 2108 1288
## m0[5] 6.78 6.75 2.58 2.51 2.49 10.89 1.00 1923 1557
## m0[6] 11.78 11.84 3.79 3.78 5.72 17.81 1.00 2001 1628
## m0[7] 22.68 22.64 4.38 4.41 15.31 29.80 1.01 2151 1382
## m0[8] 16.61 16.58 3.46 3.34 10.98 22.20 1.00 1945 1250
## m0[9] 36.75 36.83 2.16 2.06 33.06 40.23 1.00 1590 1646
## m0[10] 67.22 67.21 2.18 2.14 63.68 70.87 1.01 1862 1209
## m0[11] 24.18 24.13 2.46 2.40 20.16 28.10 1.00 1956 1358
## m0[12] 74.38 74.33 4.34 4.39 67.44 81.31 1.00 2219 1529
## m0[13] 70.49 70.45 2.37 2.36 66.63 74.41 1.01 1875 1366
## m0[14] 33.71 33.59 3.43 3.40 28.10 39.34 1.00 1799 1249
## m0[15] 24.12 24.16 1.60 1.59 21.42 26.72 1.01 2323 1377
## m0[16] 76.42 76.45 4.29 4.16 69.14 83.46 1.00 402 654
## m0[17] 29.85 29.86 2.19 2.11 26.27 33.53 1.00 2267 1784
## m0[18] 25.38 25.37 1.81 1.73 22.45 28.48 1.00 2515 1494
## m0[19] 69.54 69.50 2.24 2.17 65.97 73.26 1.00 2134 1482
## m0[20] 53.75 53.76 2.94 2.80 48.77 58.64 1.00 2228 1831
## m0[21] 11.71 11.76 3.48 3.48 5.98 17.50 1.00 1875 1438
## m0[22] 36.40 36.39 1.29 1.31 34.25 38.38 1.00 2193 1633
## m0[23] 22.49 22.53 1.66 1.67 19.70 25.20 1.00 2339 1320
## m0[24] 67.46 67.43 2.13 2.06 64.04 70.99 1.00 2240 1481
## m0[25] 15.77 15.74 2.72 2.57 11.36 20.36 1.00 2849 1555
## m0[26] 37.92 37.89 1.69 1.64 35.13 40.64 1.00 2158 1661
## m0[27] 18.09 18.06 2.40 2.31 14.17 22.06 1.00 2443 1154
## m0[28] 53.21 53.22 3.11 3.10 48.30 58.21 1.00 1729 1347
## m0[29] 40.01 40.04 1.84 1.77 36.88 43.00 1.00 1872 1657
## m0[30] -0.67 -0.61 4.01 3.95 -7.09 5.98 1.00 2514 1424
## m0[31] 29.21 29.26 1.43 1.45 26.82 31.50 1.01 2268 1649
## m0[32] 36.09 36.09 1.64 1.62 33.38 38.76 1.00 2129 1624
## m0[33] 50.51 50.43 2.92 2.97 45.79 55.21 1.00 2108 1593
## m0[34] 49.33 49.31 2.43 2.43 45.41 53.31 1.00 1888 1548
## m0[35] 20.05 20.07 2.22 2.16 16.41 23.78 1.00 2433 1171
## m0[36] 39.98 39.97 1.28 1.27 37.91 42.00 1.00 2162 1671
## m0[37] 40.13 40.09 1.80 1.74 37.14 43.06 1.00 2181 1559
## m0[38] 27.37 27.38 2.42 2.44 23.33 31.29 1.00 1830 1074
## m0[39] 13.59 13.64 2.06 2.07 10.22 16.99 1.00 2385 1199
## m0[40] 26.05 26.04 1.77 1.70 23.20 29.07 1.00 2493 1479
## m0[41] 37.58 37.56 1.63 1.58 34.89 40.31 1.00 2071 1571
## m0[42] 29.38 29.53 3.41 3.37 23.50 34.73 1.00 1322 1022
## m0[43] 70.71 70.67 2.30 2.24 67.02 74.48 1.00 2132 1547
## m0[44] 17.94 17.91 2.46 2.49 13.96 21.99 1.00 1915 1619
## m0[45] 19.00 18.97 2.55 2.43 14.93 23.26 1.00 2751 1583
## m0[46] 62.36 62.39 1.98 1.98 59.10 65.55 1.01 1892 1194
## m0[47] 90.62 90.61 4.53 4.42 83.09 97.87 1.01 570 816
## m0[48] 33.37 33.32 1.54 1.47 30.83 36.06 1.00 2226 1576
## m0[49] 64.00 64.04 2.03 2.02 60.64 67.35 1.01 1865 1194
## m0[50] 41.07 41.00 1.87 1.83 38.04 44.12 1.00 2180 1557
## y0[1] 17.19 17.12 4.74 4.59 9.51 24.99 1.00 1836 1705
## y0[2] 51.59 51.47 5.08 5.09 43.38 59.91 1.00 2035 1890
## y0[3] 53.53 53.66 4.53 4.50 46.42 60.91 1.00 2035 1972
## y0[4] 43.05 43.16 4.62 4.59 35.26 50.20 1.00 1794 1974
## y0[5] 6.65 6.52 5.18 4.84 -1.85 15.45 1.00 1831 1714
## y0[6] 11.86 12.01 5.65 5.40 2.80 20.95 1.00 1873 1374
## y0[7] 22.67 22.64 6.05 5.87 12.88 32.87 1.00 2248 1767
## y0[8] 16.58 16.69 5.59 5.57 7.45 25.68 1.00 2001 1963
## y0[9] 36.90 36.90 4.79 4.83 29.16 44.62 1.00 1779 1842
## y0[10] 67.20 67.14 4.91 4.97 59.21 75.29 1.00 1902 1951
## y0[11] 24.05 24.01 4.86 4.77 16.20 31.84 1.00 2004 1752
## y0[12] 74.36 74.32 6.13 6.10 64.45 84.38 1.00 1931 1931
## y0[13] 70.35 70.23 5.01 4.85 62.31 78.61 1.00 1991 1930
## y0[14] 33.85 33.98 5.54 5.61 24.82 42.89 1.00 1973 1957
## y0[15] 24.35 24.37 4.56 4.54 16.80 31.99 1.00 2009 1931
## y0[16] 76.46 76.41 6.19 6.09 66.71 87.00 1.00 760 1472
## y0[17] 29.82 29.82 4.80 4.65 21.88 37.81 1.00 1962 1931
## y0[18] 25.20 25.30 4.58 4.48 17.42 32.69 1.00 2036 1888
## y0[19] 69.53 69.64 4.80 4.87 61.48 76.99 1.00 1978 1760
## y0[20] 53.72 53.69 5.03 5.06 45.46 62.00 1.00 1999 1786
## y0[21] 11.62 11.57 5.59 5.12 2.45 20.92 1.00 1997 1796
## y0[22] 36.43 36.42 4.61 4.56 28.97 43.95 1.00 2083 2007
## y0[23] 22.59 22.54 4.68 4.67 15.05 30.25 1.00 2074 1861
## y0[24] 67.25 67.33 4.88 4.65 59.21 75.30 1.00 2118 2105
## y0[25] 15.76 15.77 4.97 4.95 7.79 23.86 1.00 2073 1881
## y0[26] 37.80 37.79 4.65 4.64 30.04 45.62 1.00 2079 2051
## y0[27] 18.12 18.06 4.97 4.94 9.73 26.24 1.00 2156 1873
## y0[28] 53.22 53.05 5.25 5.11 44.86 61.98 1.00 1922 1979
## y0[29] 40.10 40.08 4.61 4.57 32.28 47.39 1.00 1973 1943
## y0[30] -0.77 -0.76 5.91 5.89 -10.30 8.98 1.00 2112 1842
## y0[31] 29.13 29.23 4.60 4.54 21.48 36.72 1.00 1779 1825
## y0[32] 36.10 36.15 4.61 4.57 28.62 43.64 1.00 2013 1898
## y0[33] 50.68 50.90 5.18 5.22 42.23 58.68 1.00 2032 1926
## y0[34] 49.35 49.20 5.00 4.77 41.11 57.78 1.00 2221 1931
## y0[35] 20.25 20.20 5.04 4.84 12.09 28.56 1.00 2146 2056
## y0[36] 39.92 39.98 4.68 4.43 32.38 47.64 1.00 2083 1933
## y0[37] 40.15 40.06 4.74 4.75 32.56 48.12 1.00 2138 1816
## y0[38] 27.42 27.53 4.88 4.83 19.46 35.16 1.00 1559 1785
## y0[39] 13.49 13.56 4.91 4.71 5.39 21.43 1.00 2039 1955
## y0[40] 26.01 25.92 4.73 4.69 18.39 33.69 1.00 2104 1913
## y0[41] 37.63 37.81 4.72 4.70 29.37 45.38 1.00 2003 1695
## y0[42] 29.18 29.12 5.61 5.39 19.85 38.20 1.00 1623 1678
## y0[43] 70.80 70.82 4.85 4.81 62.66 79.01 1.00 1997 1892
## y0[44] 17.95 17.93 5.16 5.07 9.59 26.45 1.00 1997 1781
## y0[45] 18.91 18.96 5.11 5.04 10.39 27.44 1.00 2093 1767
## y0[46] 62.35 62.36 4.74 4.52 54.37 70.01 1.00 2011 1970
## y0[47] 90.65 90.71 6.28 6.26 80.12 100.80 1.01 855 1486
## y0[48] 33.27 33.32 4.67 4.65 25.54 40.77 1.00 1894 1881
## y0[49] 64.04 63.93 4.74 4.62 56.32 72.11 1.00 2022 1810
## y0[50] 40.83 40.91 4.91 4.93 32.80 48.77 1.00 1988 1846
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
estimate as class have relation
all b0l,b1l have a distribution
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)ex12-2.stan
class have relation
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
real b00;
real<lower=0,upper=100> sb0;
vector[K] b0;
real b10;
real<lower=0,upper=100> sb1;
vector[K] b1;
real<lower=0,upper=100> s;
}
model{
b0~normal(b00,sb0);
b1~normal(b10,sb1);
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-2.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -550.305
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -58.4879 0.000681448 169.74 1 1 133
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 34.9669 0.0759356 2.42048e+12 0.4956 0.4956 308
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 237 88.3164 0.0728787 1.11027e+16 1e-12 0.001 434 LS failed, Hessian reset
## Finished in 0.1 seconds.try(print(mle))## Error : Fitting failed. Unable to print.mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 1 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 1 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 2 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 2 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 3 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 3 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 4 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 4 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 2 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 1 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 3 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 4 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 2 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 2 finished in 1.3 seconds.
## Chain 1 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 3 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 4 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 1 finished in 1.5 seconds.
## Chain 3 finished in 1.4 seconds.
## Chain 4 finished in 1.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.4 seconds.
## Total execution time: 1.5 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -120.13 -119.90 4.70 4.33 -127.97 -113.43 1.01 320 77
## b00 8.47 8.42 3.21 2.76 3.37 13.55 1.00 6940 5384
## sb0 7.12 6.54 3.78 3.07 2.17 14.02 1.02 286 69
## b0[1] 0.05 0.01 4.94 5.05 -7.99 8.11 1.01 853 1565
## b0[2] 4.44 4.61 4.43 4.45 -3.24 11.26 1.00 4331 5670
## b0[3] 12.90 11.80 8.00 6.72 1.64 27.29 1.00 6395 4832
## b0[4] 6.78 6.91 3.02 2.91 1.71 11.64 1.00 7355 6071
## b0[5] 13.35 13.14 4.19 4.25 6.91 20.46 1.00 1210 5644
## b0[6] 5.54 5.78 3.89 3.79 -1.21 11.70 1.00 6341 6273
## b0[7] 11.64 11.64 2.59 2.60 7.48 15.90 1.00 1307 5860
## b0[8] 9.12 9.02 4.74 4.46 1.38 17.07 1.00 7109 5692
## b0[9] 12.38 12.43 2.18 2.15 8.72 15.95 1.00 784 883
## b10 3.02 3.02 0.48 0.45 2.25 3.78 1.00 5362 5121
## sb1 1.31 1.23 0.45 0.37 0.76 2.12 1.00 5395 4455
## b1[1] 0.77 0.77 0.43 0.44 0.05 1.48 1.01 740 794
## b1[2] 3.49 3.47 0.48 0.47 2.72 4.31 1.00 4562 5129
## b1[3] 4.32 4.37 0.51 0.46 3.39 5.07 1.00 6072 4744
## b1[4] 3.85 3.86 0.34 0.33 3.30 4.41 1.00 7084 6019
## b1[5] 3.37 3.38 0.31 0.31 2.86 3.85 1.00 1303 3245
## b1[6] 3.30 3.28 0.69 0.68 2.22 4.46 1.00 6042 5911
## b1[7] 2.66 2.65 0.28 0.28 2.20 3.11 1.00 1540 3804
## b1[8] 2.40 2.41 0.35 0.34 1.83 2.95 1.00 5427 5524
## b1[9] 2.95 2.94 0.18 0.18 2.65 3.24 1.00 967 2826
## s 4.30 4.24 0.55 0.53 3.49 5.29 1.00 5213 5021
## m0[1] 16.60 16.64 1.98 1.95 13.32 19.85 1.00 818 1484
## m0[2] 52.04 52.02 2.59 2.55 47.82 56.30 1.00 6158 4052
## m0[3] 53.45 53.47 1.56 1.48 50.88 55.96 1.00 8672 5249
## m0[4] 43.15 43.14 1.76 1.69 40.28 46.09 1.00 2092 638
## m0[5] 8.20 8.14 2.55 2.45 4.06 12.47 1.00 7573 5451
## m0[6] 12.72 12.65 3.85 3.87 6.47 19.13 1.00 2292 1526
## m0[7] 20.17 19.99 3.65 3.74 14.62 26.44 1.00 1316 5540
## m0[8] 19.15 19.23 2.63 2.63 14.66 23.29 1.00 4307 6104
## m0[9] 36.28 36.28 1.75 1.72 33.39 39.14 1.00 4419 5294
## m0[10] 67.24 67.24 2.18 2.13 63.72 70.84 1.00 7276 5694
## m0[11] 25.77 25.81 1.97 1.97 22.49 28.93 1.00 4885 6390
## m0[12] 73.23 73.29 4.20 4.11 66.21 80.05 1.00 3895 3005
## m0[13] 70.71 70.68 2.36 2.33 66.89 74.58 1.00 7049 5901
## m0[14] 32.20 32.18 3.03 2.96 27.28 37.10 1.00 7700 5950
## m0[15] 23.54 23.55 1.68 1.65 20.79 26.25 1.00 959 1928
## m0[16] 73.95 73.93 2.89 2.86 69.28 78.69 1.00 6861 5956
## m0[17] 30.27 30.28 2.12 2.09 26.78 33.73 1.00 2790 4921
## m0[18] 24.80 24.78 1.66 1.66 22.16 27.53 1.00 3415 6575
## m0[19] 69.85 69.85 2.28 2.18 66.14 73.56 1.00 5340 6021
## m0[20] 54.25 54.22 2.78 2.78 49.69 58.82 1.00 5758 4152
## m0[21] 13.81 13.87 2.59 2.52 9.45 17.96 1.00 7321 6577
## m0[22] 36.06 36.07 1.34 1.34 33.88 38.24 1.00 2492 4813
## m0[23] 21.88 21.90 1.75 1.72 19.01 24.71 1.00 917 1471
## m0[24] 67.73 67.74 2.17 2.07 64.19 71.27 1.00 5682 6109
## m0[25] 16.67 16.74 2.48 2.39 12.48 20.70 1.00 6312 5615
## m0[26] 37.77 37.76 1.65 1.60 35.06 40.46 1.00 7586 5255
## m0[27] 17.22 17.21 2.14 2.14 13.80 20.74 1.00 1580 6154
## m0[28] 53.94 53.92 2.68 2.61 49.58 58.41 1.00 2227 641
## m0[29] 39.78 39.78 1.68 1.65 36.99 42.58 1.00 2706 1312
## m0[30] 1.46 1.42 4.27 4.34 -5.49 8.37 1.01 960 1767
## m0[31] 28.73 28.73 1.50 1.49 26.31 31.16 1.00 1179 4341
## m0[32] 36.17 36.18 1.60 1.56 33.56 38.79 1.00 7522 5110
## m0[33] 48.78 48.77 2.56 2.49 44.65 52.99 1.00 6446 4969
## m0[34] 49.78 49.75 2.24 2.16 46.13 53.51 1.00 2006 633
## m0[35] 19.26 19.25 1.99 2.00 16.07 22.54 1.00 1821 6554
## m0[36] 39.70 39.69 1.32 1.30 37.54 41.82 1.00 4850 4724
## m0[37] 39.70 39.69 1.75 1.70 36.82 42.57 1.00 7524 5419
## m0[38] 26.90 26.88 2.34 2.29 23.11 30.75 1.00 8121 5468
## m0[39] 12.80 12.85 2.16 2.14 9.18 16.33 1.00 784 866
## m0[40] 25.49 25.47 1.63 1.61 22.91 28.17 1.00 3954 6891
## m0[41] 37.47 37.46 1.60 1.52 34.86 40.11 1.00 9189 5003
## m0[42] 28.38 28.35 2.35 2.26 24.58 32.18 1.00 7522 6148
## m0[43] 71.05 71.04 2.34 2.25 67.23 74.86 1.00 5171 6182
## m0[44] 19.01 19.08 2.10 2.05 15.54 22.50 1.00 7791 5660
## m0[45] 19.79 19.84 2.35 2.29 15.84 23.62 1.00 5408 6091
## m0[46] 62.11 62.12 1.99 1.94 58.92 65.39 1.00 7457 5149
## m0[47] 91.78 91.79 3.47 3.46 85.99 97.46 1.00 6769 6210
## m0[48] 33.10 33.07 1.49 1.42 30.66 35.56 1.00 9143 5179
## m0[49] 63.85 63.86 2.04 2.00 60.57 67.22 1.00 7425 5261
## m0[50] 40.52 40.51 1.80 1.75 37.55 43.49 1.00 7474 5068
## y0[1] 16.71 16.74 4.75 4.62 8.77 24.52 1.00 6739 7967
## y0[2] 52.07 52.04 5.02 4.85 43.87 60.32 1.00 7624 7790
## y0[3] 53.43 53.44 4.61 4.53 45.90 60.91 1.00 8458 7669
## y0[4] 43.15 43.17 4.60 4.62 35.63 50.70 1.00 6986 7610
## y0[5] 8.17 8.12 5.02 4.99 0.01 16.46 1.00 7784 7191
## y0[6] 12.69 12.70 5.76 5.68 3.22 22.02 1.00 6821 7247
## y0[7] 20.17 20.11 5.67 5.61 10.93 29.43 1.00 4865 6320
## y0[8] 19.19 19.21 5.07 4.96 10.70 27.59 1.00 6405 6878
## y0[9] 36.29 36.32 4.67 4.47 28.51 43.90 1.00 8123 7924
## y0[10] 67.26 67.28 4.78 4.63 59.26 75.14 1.00 7549 7593
## y0[11] 25.79 25.76 4.76 4.67 18.08 33.58 1.00 7490 7807
## y0[12] 73.25 73.35 6.07 5.88 63.44 82.92 1.00 6285 6965
## y0[13] 70.68 70.71 4.94 4.83 62.67 78.80 1.00 7997 7496
## y0[14] 32.22 32.17 5.25 5.13 23.52 40.78 1.00 8116 7492
## y0[15] 23.53 23.57 4.63 4.55 15.94 30.95 1.00 7439 7574
## y0[16] 73.92 73.96 5.28 5.22 65.27 82.64 1.00 7296 6870
## y0[17] 30.30 30.21 4.83 4.74 22.44 38.11 1.00 7869 7133
## y0[18] 24.88 24.88 4.70 4.64 17.23 32.47 1.00 8104 7953
## y0[19] 69.88 69.89 4.90 4.79 61.82 77.94 1.00 8002 6998
## y0[20] 54.25 54.26 5.12 5.13 46.03 62.82 1.00 7878 7733
## y0[21] 13.79 13.78 5.04 4.98 5.55 22.20 1.00 7359 7637
## y0[22] 36.01 36.05 4.59 4.48 28.42 43.43 1.00 7823 7522
## y0[23] 21.80 21.83 4.64 4.54 14.24 29.38 1.00 5683 7161
## y0[24] 67.70 67.68 4.91 4.78 59.74 75.84 1.00 8215 7063
## y0[25] 16.72 16.71 4.99 4.93 8.62 24.87 1.00 8032 7907
## y0[26] 37.77 37.75 4.62 4.52 30.31 45.46 1.00 8221 7875
## y0[27] 17.30 17.27 4.86 4.68 9.44 25.41 1.00 7736 7325
## y0[28] 53.93 53.93 5.08 4.92 45.61 62.31 1.00 6234 5673
## y0[29] 39.71 39.70 4.63 4.53 32.02 47.32 1.00 7320 7737
## y0[30] 1.49 1.41 6.05 6.02 -8.30 11.48 1.00 3045 5252
## y0[31] 28.83 28.84 4.60 4.42 21.37 36.43 1.00 7363 7283
## y0[32] 36.17 36.13 4.59 4.53 28.69 43.70 1.00 7819 7187
## y0[33] 48.75 48.85 5.07 4.99 40.34 57.01 1.00 8080 7857
## y0[34] 49.71 49.71 4.89 4.85 41.67 57.59 1.00 6906 7733
## y0[35] 19.27 19.36 4.87 4.72 11.24 27.07 1.00 6420 7923
## y0[36] 39.72 39.72 4.60 4.54 32.15 47.31 1.00 7721 7491
## y0[37] 39.70 39.72 4.63 4.54 32.09 47.19 1.00 7999 7803
## y0[38] 26.91 26.91 4.95 4.86 19.00 35.14 1.00 8000 7650
## y0[39] 12.86 12.88 4.89 4.73 4.81 20.93 1.00 6750 7059
## y0[40] 25.45 25.48 4.70 4.62 17.66 33.11 1.00 7943 8014
## y0[41] 37.45 37.43 4.63 4.59 29.94 44.99 1.00 8094 7680
## y0[42] 28.44 28.40 4.85 4.75 20.55 36.34 1.00 8191 7636
## y0[43] 71.07 71.06 4.95 4.96 63.12 79.23 1.00 8041 7913
## y0[44] 19.04 19.00 4.76 4.66 11.23 26.92 1.00 8238 7365
## y0[45] 19.80 19.83 4.95 4.92 11.69 27.96 1.00 7790 7353
## y0[46] 62.15 62.12 4.81 4.70 54.27 70.02 1.00 7777 7664
## y0[47] 91.75 91.74 5.60 5.58 82.48 100.86 1.00 7896 7053
## y0[48] 33.02 33.04 4.57 4.47 25.54 40.55 1.00 7887 7521
## y0[49] 63.82 63.89 4.81 4.67 55.84 71.64 1.00 7305 7229
## y0[50] 40.47 40.44 4.70 4.60 32.77 48.19 1.00 8113 7939
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
generalized linear mixed model
(X,y)i=1-n
b[b0,b1,...]
linear model y~N(Xb,s)
generalized linear model y~dist.(m=link(Xb),s)
fixed effect b0, b1,...
individual random effect b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...
class c=1-k
class effect b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...
for y=dist.(m,s)
random intercept model m=b0+r0i+b1*x, m=b0+r0c+b1*x
random coefficient model m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x
mixed model m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x
note
@ yi=b0+b1*xi+r0i is not useful, ri is included in s
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
but variance is larger than mean (over dispersion),
random intercept model is usefuln=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)
ex13-0.stan
generalized linear model, poisson regression
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~poisson_log(b0+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-0.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -4098.1
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 1479.7 0.000581323 0.00909063 1 1 19
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 1479.70
## b0 2.09
## b1 0.21
## m0[1] 2.31
## m0[2] 2.89
## m0[3] 3.45
## m0[4] 2.70
## m0[5] 3.41
## m0[6] 2.71
## m0[7] 3.67
## m0[8] 3.78
## m0[9] 2.89
## m0[10] 3.64
## m0[11] 3.77
## m0[12] 3.85
## m0[13] 3.78
## m0[14] 3.77
## m0[15] 3.13
## m0[16] 3.48
## m0[17] 2.67
## m0[18] 3.08
## m0[19] 2.46
## m0[20] 3.97
## y0[1] 5.00
## y0[2] 13.00
## y0[3] 55.00
## y0[4] 22.00
## y0[5] 34.00
## y0[6] 22.00
## y0[7] 48.00
## y0[8] 42.00
## y0[9] 17.00
## y0[10] 29.00
## y0[11] 45.00
## y0[12] 52.00
## y0[13] 43.00
## y0[14] 42.00
## y0[15] 18.00
## y0[16] 35.00
## y0[17] 19.00
## y0[18] 17.00
## y0[19] 4.00
## y0[20] 59.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 1478.72 1479.05 1.00 0.68 1476.68 1479.65 1.01 598 762
## b0 2.10 2.10 0.14 0.13 1.88 2.35 1.02 228 209
## b1 0.21 0.21 0.02 0.02 0.17 0.24 1.02 236 249
## m0[1] 2.32 2.32 0.12 0.11 2.13 2.53 1.02 230 226
## m0[2] 2.89 2.89 0.07 0.07 2.78 3.01 1.01 266 291
## m0[3] 3.45 3.46 0.04 0.04 3.39 3.52 1.00 1271 1359
## m0[4] 2.71 2.71 0.08 0.08 2.58 2.86 1.02 246 220
## m0[5] 3.41 3.41 0.04 0.04 3.34 3.48 1.00 971 1379
## m0[6] 2.72 2.72 0.08 0.08 2.58 2.86 1.02 246 220
## m0[7] 3.67 3.67 0.04 0.04 3.60 3.75 1.00 1755 1523
## m0[8] 3.78 3.78 0.05 0.05 3.69 3.86 1.00 1219 1598
## m0[9] 2.90 2.90 0.07 0.07 2.79 3.02 1.01 267 291
## m0[10] 3.64 3.64 0.04 0.04 3.57 3.71 1.00 1890 1533
## m0[11] 3.77 3.77 0.05 0.05 3.69 3.85 1.00 1259 1598
## m0[12] 3.85 3.85 0.05 0.05 3.76 3.93 1.00 866 1185
## m0[13] 3.78 3.78 0.05 0.05 3.70 3.86 1.00 1176 1546
## m0[14] 3.77 3.77 0.05 0.05 3.68 3.85 1.00 1275 1602
## m0[15] 3.13 3.13 0.05 0.05 3.05 3.22 1.01 348 498
## m0[16] 3.49 3.49 0.04 0.04 3.42 3.55 1.00 1447 1360
## m0[17] 2.68 2.68 0.09 0.08 2.54 2.83 1.02 244 231
## m0[18] 3.08 3.08 0.06 0.05 2.99 3.18 1.01 319 367
## m0[19] 2.47 2.47 0.10 0.10 2.31 2.66 1.02 234 239
## m0[20] 3.97 3.97 0.06 0.06 3.87 4.07 1.01 585 1000
## y0[1] 10.28 10.00 3.52 2.97 5.00 16.00 1.00 1071 1822
## y0[2] 18.13 18.00 4.47 4.45 11.00 26.00 1.00 1632 1856
## y0[3] 31.50 31.00 5.83 5.93 22.00 41.00 1.00 2029 1947
## y0[4] 15.26 15.00 4.10 4.45 9.00 22.00 1.00 1506 1754
## y0[5] 30.32 30.00 5.57 5.93 22.00 40.00 1.00 1545 1886
## y0[6] 15.09 15.00 4.15 4.45 9.00 22.00 1.00 1176 1574
## y0[7] 39.56 39.00 6.76 7.41 29.00 51.00 1.00 2109 2105
## y0[8] 43.67 43.00 7.09 7.41 32.00 56.00 1.00 1739 1715
## y0[9] 18.21 18.00 4.34 4.45 11.00 26.00 1.00 1724 1839
## y0[10] 38.11 38.00 6.09 5.93 29.00 48.00 1.00 2009 1629
## y0[11] 43.51 43.00 6.93 7.41 33.00 55.00 1.00 2154 2016
## y0[12] 46.92 47.00 7.44 7.41 35.00 60.00 1.00 1979 1889
## y0[13] 43.90 44.00 6.98 7.41 33.00 56.00 1.00 1948 2023
## y0[14] 43.25 43.00 7.11 7.41 32.00 55.00 1.00 1826 1931
## y0[15] 22.94 23.00 4.81 4.45 15.00 31.00 1.00 1659 1822
## y0[16] 32.55 32.00 5.96 5.93 23.00 43.00 1.00 2112 1810
## y0[17] 14.53 14.00 3.96 4.45 8.00 22.00 1.00 1342 1960
## y0[18] 21.74 22.00 4.75 4.45 14.00 30.00 1.00 1684 1927
## y0[19] 12.00 12.00 3.70 2.97 6.00 18.00 1.00 1549 1894
## y0[20] 52.72 53.00 7.91 7.41 40.00 66.00 1.00 1364 1839
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
ex13-1.stan
generalized linear mixed model, poisson regression
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
real<lower=0> sr0;
vector[N] r0;
}
model{
r0~normal(0,sr0);
for(i in 1:N)
y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+r0[i]+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-1.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = 4594.93
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 29722.6 0.149664 4711.66 0.4869 0.7899 155
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 29874.8 0.193417 3.59653e+12 0.4627 0.4627 325
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 213 29877.3 0.00213885 2.78206e+11 0.6069 0.6069 343
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 29877.30
## b0 1.38
## b1 0.32
## sr0 0.00
## r0[1] 0.00
## r0[2] 0.00
## r0[3] 0.00
## r0[4] 0.00
## r0[5] 0.00
## r0[6] 0.00
## r0[7] 0.00
## r0[8] 0.00
## r0[9] 0.00
## r0[10] 0.00
## r0[11] 0.00
## r0[12] 0.00
## r0[13] 0.00
## r0[14] 0.00
## r0[15] 0.00
## r0[16] 0.00
## r0[17] 0.00
## r0[18] 0.00
## r0[19] 0.00
## r0[20] 0.00
## m0[1] 1.71
## m0[2] 2.61
## m0[3] 3.49
## m0[4] 2.33
## m0[5] 3.42
## m0[6] 2.33
## m0[7] 3.83
## m0[8] 3.99
## m0[9] 2.62
## m0[10] 3.78
## m0[11] 3.99
## m0[12] 4.11
## m0[13] 4.00
## m0[14] 3.98
## m0[15] 2.99
## m0[16] 3.54
## m0[17] 2.28
## m0[18] 2.90
## m0[19] 1.96
## m0[20] 4.30
## y0[1] 7.00
## y0[2] 14.00
## y0[3] 35.00
## y0[4] 13.00
## y0[5] 32.00
## y0[6] 13.00
## y0[7] 50.00
## y0[8] 56.00
## y0[9] 18.00
## y0[10] 40.00
## y0[11] 65.00
## y0[12] 51.00
## y0[13] 46.00
## y0[14] 45.00
## y0[15] 11.00
## y0[16] 37.00
## y0[17] 13.00
## y0[18] 11.00
## y0[19] 5.00
## y0[20] 84.00mcmc=goMCMC(mdl,data,wrm=500)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2 finished in 1.3 seconds.
## Chain 3 finished in 1.2 seconds.
## Chain 4 finished in 1.2 seconds.
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 1.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.6 seconds.seeMCMC(mcmc,'[m,r]')## variable mean median sd mad q5 q95 rhat ess_bulk
## lp__ 29683.39 29682.80 15.22 16.01 29658.60 29709.90 1.11 38
## b0 2.09 2.09 0.03 0.03 2.04 2.15 1.01 307
## b1 0.21 0.21 0.00 0.00 0.20 0.22 1.01 324
## sr0 0.01 0.01 0.01 0.00 0.00 0.02 1.11 34
## r0[1] 0.00 0.00 0.01 0.00 -0.01 0.01 1.03 2958
## r0[2] 0.00 0.00 0.01 0.00 -0.01 0.01 1.05 3644
## r0[3] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3359
## r0[4] 0.00 0.00 0.01 0.00 -0.01 0.01 1.05 3047
## r0[5] 0.00 0.00 0.01 0.00 -0.01 0.01 1.05 3155
## r0[6] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3246
## r0[7] 0.00 0.00 0.01 0.00 -0.01 0.01 1.03 3608
## r0[8] 0.00 0.00 0.01 0.00 -0.01 0.01 1.05 3414
## r0[9] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 2472
## r0[10] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3250
## r0[11] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3013
## r0[12] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3278
## r0[13] 0.00 0.00 0.01 0.00 -0.01 0.01 1.02 3034
## r0[14] 0.00 0.00 0.01 0.00 -0.01 0.01 1.03 2981
## r0[15] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 2324
## r0[16] 0.00 0.00 0.01 0.00 -0.01 0.01 1.05 3144
## r0[17] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3057
## r0[18] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3203
## r0[19] 0.00 0.00 0.01 0.00 -0.01 0.01 1.04 3609
## r0[20] 0.00 0.00 0.01 0.00 -0.01 0.01 1.05 3338
## m0[1] 2.31 2.31 0.03 0.03 2.26 2.36 1.01 344
## m0[2] 2.89 2.89 0.02 0.02 2.86 2.92 1.01 408
## m0[3] 3.45 3.45 0.01 0.01 3.43 3.47 1.00 1879
## m0[4] 2.70 2.70 0.02 0.02 2.67 2.74 1.01 362
## m0[5] 3.41 3.41 0.01 0.01 3.39 3.43 1.00 1555
## m0[6] 2.71 2.71 0.02 0.02 2.67 2.75 1.00 365
## m0[7] 3.67 3.67 0.01 0.01 3.65 3.69 1.00 2375
## m0[8] 3.78 3.78 0.01 0.01 3.75 3.80 1.00 1994
## m0[9] 2.89 2.89 0.02 0.02 2.86 2.92 1.01 407
## m0[10] 3.64 3.64 0.01 0.01 3.62 3.66 1.00 2545
## m0[11] 3.77 3.77 0.01 0.01 3.75 3.79 1.01 1670
## m0[12] 3.85 3.85 0.01 0.01 3.83 3.87 1.00 1390
## m0[13] 3.78 3.78 0.01 0.01 3.76 3.81 1.00 1787
## m0[14] 3.77 3.77 0.01 0.01 3.75 3.79 1.00 1660
## m0[15] 3.13 3.13 0.01 0.01 3.11 3.15 1.00 545
## m0[16] 3.48 3.48 0.01 0.01 3.47 3.50 1.00 2215
## m0[17] 2.67 2.67 0.02 0.02 2.64 2.71 1.01 372
## m0[18] 3.08 3.08 0.02 0.02 3.05 3.10 1.00 542
## m0[19] 2.47 2.47 0.03 0.03 2.42 2.51 1.00 332
## m0[20] 3.97 3.97 0.02 0.02 3.95 4.00 1.00 928
## y0[1] 10.09 10.00 3.23 2.97 5.00 16.00 1.00 1957
## y0[2] 17.96 18.00 4.34 4.45 11.00 25.00 1.00 1712
## y0[3] 31.46 31.00 5.79 5.93 22.00 41.05 1.00 2015
## y0[4] 14.91 15.00 4.00 4.45 9.00 22.00 1.00 1896
## y0[5] 30.26 30.00 5.54 5.93 21.00 39.05 1.00 1897
## y0[6] 15.03 15.00 3.83 4.45 9.00 21.05 1.00 1874
## y0[7] 39.36 39.00 6.12 5.93 30.00 50.00 1.00 1970
## y0[8] 43.48 43.00 6.58 5.93 33.00 55.00 1.00 2050
## y0[9] 17.92 18.00 4.19 4.45 11.00 25.00 1.00 1884
## y0[10] 38.13 38.00 6.48 5.93 28.00 49.00 1.00 1990
## y0[11] 43.38 43.00 6.68 7.41 32.00 55.00 1.00 1743
## y0[12] 46.89 47.00 6.94 7.41 36.00 59.00 1.00 1918
## y0[13] 43.97 44.00 6.65 5.93 33.00 55.00 1.00 1974
## y0[14] 43.39 43.00 6.58 5.93 33.00 54.00 1.00 2138
## y0[15] 23.09 23.00 4.85 4.45 15.00 31.00 1.00 2040
## y0[16] 32.45 32.00 5.71 5.93 23.00 42.00 1.00 2013
## y0[17] 14.56 14.00 3.80 4.45 8.00 21.00 1.00 1949
## y0[18] 21.53 21.00 4.71 4.45 14.00 30.00 1.00 2139
## y0[19] 11.71 11.00 3.38 2.97 7.00 17.05 1.00 1909
## y0[20] 53.23 53.00 7.26 7.41 42.00 65.00 1.00 1886
## ess_tail
## 78
## 377
## 329
## 56
## 581
## 598
## 687
## 765
## 600
## 906
## 527
## 478
## 517
## 605
## 735
## 554
## 527
## 771
## 580
## 824
## 557
## 663
## 656
## 678
## 474
## 686
## 755
## 587
## 1012
## 607
## 1018
## 793
## 993
## 980
## 1200
## 1139
## 1278
## 1386
## 1080
## 1039
## 583
## 1074
## 477
## 833
## 1941
## 1704
## 1675
## 1771
## 1936
## 1774
## 1953
## 1906
## 1895
## 1994
## 1799
## 1905
## 2025
## 2026
## 1985
## 1882
## 2050
## 1912
## 1858
## 1912
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
categorical dist. with parameters following Dirichlet dist.
h~Dir(a), h[0,1], sum(h)=1, a[0,1], sum(a)=1
y~Cat(h), y[0,1], sum(y)=1ex15-1.stan
data {
int N;
int K;
array[N] int<lower=1,upper=K> y;
}
parameters {
simplex[K] a;
simplex[K] h;
}
model {
h~dirichlet(a);
y~categorical(h);
}library(gtools)
n=20
h=rdirichlet(n,c(0.5,0.3,0.2))
summary(h)## V1 V2 V3
## Min. :0.001 Min. :0.000 Min. :0.000
## 1st Qu.:0.148 1st Qu.:0.032 1st Qu.:0.003
## Median :0.425 Median :0.082 Median :0.080
## Mean :0.469 Mean :0.295 Mean :0.235
## 3rd Qu.:0.791 3rd Qu.:0.696 3rd Qu.:0.313
## Max. :0.996 Max. :0.999 Max. :0.927c0=c(1,2,3)
for(i in 1:n) y[i]=sample(c0,1,prob=h[i,])
table(y) |> prop.table()## y
## 1 2 3
## 0.35 0.35 0.30data=list(N=n,K=ncol(h),y=y)
mdl=cmdstan_model('./ex15-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -26.0962
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -22.6736 0.000289964 3.89817e-05 0.9996 0.9996 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -22.67
## a[1] 0.34
## a[2] 0.34
## a[3] 0.32
## h[1] 0.35
## h[2] 0.35
## h[3] 0.30mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -31.41 -31.07 1.48 1.31 -34.29 -29.66 1.00 915 1217
## a[1] 0.34 0.32 0.18 0.19 0.08 0.66 1.00 1171 883
## a[2] 0.34 0.32 0.18 0.19 0.08 0.66 1.00 1130 960
## a[3] 0.32 0.30 0.17 0.18 0.08 0.63 1.00 1167 1244
## h[1] 0.35 0.34 0.10 0.10 0.19 0.53 1.00 2774 1644
## h[2] 0.35 0.35 0.10 0.11 0.19 0.53 1.00 2285 1639
## h[3] 0.30 0.29 0.10 0.10 0.15 0.47 1.00 2261 1375
beta binomial dist.
y~betaB(n,a,b): y~B(n,p), p~beta(a,b)ex15-2.stan
data {
int N;
int n;
array[N] int y;
}
parameters {
real<lower=0> a;
real<lower=0> b;
}
model {
a~cauchy(0,5);
b~cauchy(0,5);
y~beta_binomial(n,a,b);
}
generated quantities{
int y1;
y1=beta_binomial_rng(n,a,b);
}n=30
a=3
b=4
p=rbeta(n,a,b)
n0=10
y=rbinom(n,n0,p)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.25 5.00 4.43 6.00 8.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex15-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -246.925
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 -197.026 0.00201419 0.00162191 1 1 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -197.03
## a 2.10
## b 2.69
## y1 3.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -196.14 -195.76 1.19 0.77 -198.33 -195.08 1.00 577 465
## a 3.21 2.79 1.90 1.14 1.47 6.02 1.01 369 284
## b 4.09 3.58 2.42 1.44 1.89 7.38 1.01 378 304
## y1 4.39 4.00 2.38 2.97 1.00 8.00 1.00 1886 1943
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
bimodal distribution, mixed normal distribution
n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()
EM algorithm
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -239 100 6 -506 -507
##
## Clustering table:
## 1 2 3
## 32 35 33rst$parameters## $pro
## [1] 0.319 0.352 0.330
##
## $mean
## 1 2 3
## -4.7988 -0.0713 5.0721
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 0.796
##
##
## $Vinv
## NULLplot(rst)
ex17-1.stan
data {
int N;
vector[N] y;;
}
parameters {
simplex[3] h; //ratio of mix
real m1;
real m2;
real m3;
real<lower=0> s1;
real<lower=0> s2;
real<lower=0> s3;
}
model {
s1~cauchy(0,5);
s2~cauchy(0,5);
s3~cauchy(0,5);
for (i in 1:N) {
vector[3] p;
p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
target+=log_sum_exp(p);
}
}mdl=cmdstan_model('./ex17-1.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -375.528
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -275.551 0.000193854 0.0354569 1 1 126
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 102 -275.551 0.000199822 0.0235199 1 1 129
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -275.55
## h[1] 0.07
## h[2] 0.89
## h[3] 0.04
## m1 0.35
## m2 0.13
## m3 -0.52
## s1 0.17
## s2 4.30
## s3 0.02mcmc=goMCMC(mdl,data,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 3 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 4 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 3 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 2 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 4 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 4 finished in 1.4 seconds.
## Chain 1 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 3 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 3 finished in 1.6 seconds.
## Chain 2 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 2 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 2 finished in 2.1 seconds.
## Chain 1 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 1 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 1 finished in 2.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.9 seconds.
## Total execution time: 2.5 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -246.99 -246.55 2.85 2.03 -251.08 -243.98 1.00 2453 1849
## h[1] 0.34 0.33 0.05 0.05 0.26 0.42 1.03 85 1919
## h[2] 0.33 0.33 0.05 0.05 0.25 0.41 1.03 101 579
## h[3] 0.33 0.33 0.05 0.05 0.25 0.41 1.03 102 3823
## m1 0.11 -0.05 3.99 2.61 -4.92 5.21 2.41 4 26
## m2 1.34 3.55 4.10 2.86 -4.94 5.28 2.10 5 29
## m3 -1.13 -2.04 4.07 4.06 -5.03 5.21 2.33 5 27
## s1 0.99 0.91 1.48 0.16 0.71 1.27 1.16 16 89
## s2 0.95 0.92 0.21 0.15 0.73 1.25 1.15 17 71
## s3 0.99 0.96 0.19 0.17 0.73 1.32 1.22 12 76
